MODULE ADVECVAR
!
! Verwendungszweck:
!
!      Bereitstellung der Advektionsvariablen und deren Ableitungen
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   rho_adv(i) : Advektierter Werte an der Zellgrenze i
!                fuer skalare groessen:  rho   - Dichte
!                                        e     - innere Energie
!                                        J     - Strahlungsenergiedichte
!                                  bezieht sich der advektierte wert mit dem index i auf den radiuspunkt r_(i)
!                fuer vektor-groessen:   momr  - Impuls r-Komponente (rho * u_r)
!                                        momp  - Impuls phi-Komponente (rho * u_phi * r)
!                                        H     - Strahlungsfluss
!                                  bezieht sich der advektierte wert mit dem index i auf den gemittelten radiuspunkt r_(i+1/2)
!
!                demgemaess unterschieden sich auch die indizes nach denen Abgeleitet wird
!                fuer skalare groessen: (i-2), (i-1), (i), (i+1)
!                fuer vektor-groessen:         (i-1), (i), (i+1), (i+2)
!
!   drho_adv_dim2(i) : Ableitung von rho_adv nach der Dichte am Punkt i-2
!                      die anderen Ableitungen sinngemaess - auch fuer andere Advektionsgroessen
!                      Die Ableitungen werden dabei immer gegen das Argument des Advektionsschemas berechnet, d.h. XZ !
!                      Da in den Argumenten meistens Ausdrucke stecken, z.B. (rhoZ eZ), geht das auch gar nicht anders
!                      da Mathematica die innere Ableitung des Arguments berechnen koennen muss; also z.B. d(rhoZ eZ)/d(rho)
!
!-----------------------------------------------------------------------

      use config, only : rkind, np
      use physco, only : z0


      implicit none

      private
      public :: ADVEC, &
                rho_adv, drho_adv_dim2, drho_adv_dim1, drho_adv_di, drho_adv_dip1, &
                e_adv, de_adv_dim2, de_adv_dim1, de_adv_di, de_adv_dip1, &
                momp_adv, dmomp_adv_dim2, dmomp_adv_dim1, dmomp_adv_di, dmomp_adv_dip1, &
                J_adv, dJ_adv_dim2, dJ_adv_dim1, dJ_adv_di, dJ_adv_dip1, &
                momr_adv, dmomr_adv_dim1, dmomr_adv_di, dmomr_adv_dip1, dmomr_adv_dip2, &
                H_adv, dH_adv_dim1, dH_adv_di, dH_adv_dip1, dH_adv_dip2


      real(kind=rkind) :: rho_adv(np)  =z0, drho_adv_dim2(np)  =z0, drho_adv_dim1(np)  =z0, drho_adv_di(np)  =z0, drho_adv_dip1(np)  =z0
      real(kind=rkind) :: e_adv(np)    =z0, de_adv_dim2(np)    =z0, de_adv_dim1(np)    =z0, de_adv_di(np)    =z0, de_adv_dip1(np)    =z0
      real(kind=rkind) :: momp_adv(np) =z0, dmomp_adv_dim2(np) =z0, dmomp_adv_dim1(np) =z0, dmomp_adv_di(np) =z0, dmomp_adv_dip1(np) =z0
      real(kind=rkind) :: J_adv(np)    =z0, dJ_adv_dim2(np)    =z0, dJ_adv_dim1(np)    =z0, dJ_adv_di(np)    =z0, dJ_adv_dip1(np)    =z0

      real(kind=rkind) :: momr_adv(np) =z0, dmomr_adv_dim1(np) =z0, dmomr_adv_di(np) =z0, dmomr_adv_dip1(np) =z0, dmomr_adv_dip2(np) =z0
      real(kind=rkind) :: H_adv(np)    =z0, dH_adv_dim1(np)    =z0, dH_adv_di(np)    =z0, dH_adv_dip1(np)    =z0, dH_adv_dip2(np)    =z0

CONTAINS





SUBROUTINE ADVEC
!
! Verwendungszweck:
!
!   Zusammensetzung der advektierten Groessen:
!
!
!   Xout:  advektierte Groesse bei einer Ausstroemung (d.h. Stroemung nach rechts)
!   Xin :  advektierte Groesse bei einer Einstroemung (d.h. Stroemung nach links)
!
!                       |
!                     <-+-Xin
!                       |
!                  Xout-+->
!                       |
!        |         |    |    |         |
! Core --+---------+---------+---------+---> surface
!        X1        X2        X3        X4
!
!
!  Advektionsschema:   X_out = f(X1,X2,X3)
!                      X_in  = f(X2,X3,X4)
!
!  Die Groessen Xout und Xin (zusammen mit den Ableitungen) fuer EINE Zellgrenze
!  werden durch die Subroutine "adv_scheme" berechnet. Da jede Zelle zwei Zell-
!  grenzen besitzt (Oben & Unten) wird "adv_scheme" fuer jede Zelle ebenfalls
!  zweimal aufgerufen.
!
!  Es werden jeweils die 4 benachbarten Werte der zu advektierenden Groesse uebergeben.
!  d.h. also fuer die untere Grenze:
!
!     X_im2   X_im1   X_i   X_ip1  X_ip2
!      X1      X2     X3     X4
!
!  und fuer die obere Grenze:
!
!     X_im2   X_im1   X_i   X_ip1  X_ip2
!              X1      X2     X3     X4
!
!  Die Groessen Xout und Xin werden anschliessend - abhaengig von der Flussrichtung -
!  in der vorliegenden Subroutine zusammengesetzt.
!  also beispielsweise:
!
!      X_adv_O  = S_fO * Xout + (z1-S_fO) * Xin
!
!  Die Schalter fuer die Flussrichtungen S_fO ( bzw. fuer vekorielle Zellen V_fO) sind
!  jeweils bei einer Stroemung nach Aussen 1, sonst 0.
!  Die Ableitungen der advektierten Groessen werden analog zusammengestellt.
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   idx        : Schleifenindex
!   Xout       : Advektionsgroesse fuer Stroemungen nach AUSSEN
!   dXoutdX1   : Ableitung von Xout nach X1 (detto fuer X2 & X3 jedoch NICHT X4)
!   Xin        : Advektionsgroesse fuer Stroemungen nach INNEN
!   dXindX2    : Ableitung von Xin nach X2 (detto fuer X3 & X4 jedoch NICHT X1)
!
!-----------------------------------------------------------------------

      use config,  only : rkind, np
      use physco,  only : z1, z12, z0
      use geomvar, only : S_flag, V_flag
      use primvar, only : XZ, MR, MD, ME, MUr, MUp, MJ, MH


      implicit none

      integer          :: idx

      real(kind=rkind) :: Xout, dXoutdX1, dXoutdX2, dXoutdX3
      real(kind=rkind) :: Xin, dXindX2, dXindX3, dXindX4

      external ADV_SCHEME


      !$OMP PARALLEL DO PRIVATE (idx, Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)
      ! Fuer skalare Zellen brauchen wir Advektion von S_flux(3) bis S_flux(np-2)
      do idx = 3,np-2

!-----------------------------------------------------------------------
!  Advektion der Dichte
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(MD,idx-2), XZ(MD,idx-1), XZ(MD,idx), XZ(MD,idx+1), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         rho_adv(idx)        = S_flag(idx) * Xout     + (z1-S_flag(idx)) * Xin
         drho_adv_dim2(idx)  = S_flag(idx) * dXoutdX1
         drho_adv_dim1(idx)  = S_flag(idx) * dXoutdX2 + (z1-S_flag(idx)) * dXindX2
         drho_adv_di(idx)    = S_flag(idx) * dXoutdX3 + (z1-S_flag(idx)) * dXindX3
         drho_adv_dip1(idx)  =                          (z1-S_flag(idx)) * dXindX4


!-----------------------------------------------------------------------
!  Advektion der inneren Energie
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(ME,idx-2)*XZ(MD,idx-2), XZ(ME,idx-1)*XZ(MD,idx-1), XZ(ME,idx)*XZ(MD,idx), XZ(ME,idx+1)*XZ(MD,idx+1), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         e_adv(idx)          = S_flag(idx) * Xout     + (z1-S_flag(idx)) * Xin
         de_adv_dim2(idx)    = S_flag(idx) * dXoutdX1
         de_adv_dim1(idx)    = S_flag(idx) * dXoutdX2 + (z1-S_flag(idx)) * dXindX2
         de_adv_di(idx)      = S_flag(idx) * dXoutdX3 + (z1-S_flag(idx)) * dXindX3
         de_adv_dip1(idx)    =                          (z1-S_flag(idx)) * dXindX4


!-----------------------------------------------------------------------
!  Advektion des Impulses phi-Komponente
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(MD,idx-2)*XZ(MUp,idx-2), XZ(MD,idx-1)*XZ(MUp,idx-1), XZ(MD,idx)*XZ(MUp,idx), XZ(MD,idx+1)*XZ(MUp,idx+1), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         momp_adv(idx)       = S_flag(idx) * Xout     + (z1-S_flag(idx)) * Xin
         dmomp_adv_dim2(idx) = S_flag(idx) * dXoutdX1
         dmomp_adv_dim1(idx) = S_flag(idx) * dXoutdX2 + (z1-S_flag(idx)) * dXindX2
         dmomp_adv_di(idx)   = S_flag(idx) * dXoutdX3 + (z1-S_flag(idx)) * dXindX3
         dmomp_adv_dip1(idx) =                          (z1-S_flag(idx)) * dXindX4


!-----------------------------------------------------------------------
!  Advektion der Strahlungsenergiedichte
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(MJ,idx-2), XZ(MJ,idx-1), XZ(MJ,idx), XZ(MJ,idx+1), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         J_adv(idx)          = S_flag(idx) * Xout     + (z1-S_flag(idx)) * Xin
         dJ_adv_dim2(idx)    = S_flag(idx) * dXoutdX1
         dJ_adv_dim1(idx)    = S_flag(idx) * dXoutdX2 + (z1-S_flag(idx)) * dXindX2
         dJ_adv_di(idx)      = S_flag(idx) * dXoutdX3 + (z1-S_flag(idx)) * dXindX3
         dJ_adv_dip1(idx)    =                          (z1-S_flag(idx)) * dXindX4

      end do
      !$OMP END PARALLEL DO


      !$OMP PARALLEL DO PRIVATE (idx, Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)
      ! Fuer vektorielle Zellen brauchen wir Advektion von V_flux(2) bis V_flux(np-2)
      do idx = 2,np-2

!-----------------------------------------------------------------------
!  Advektion des Impulses r-Komponente
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(MD,idx-1)*XZ(MUr,idx-1), XZ(MD,idx)*XZ(MUr,idx), XZ(MD,idx+1)*XZ(MUr,idx+1), XZ(MD,idx+2)*XZ(MUr,idx+2), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         momr_adv(idx)       = V_flag(idx) * Xout     + (z1-V_flag(idx)) * Xin
         dmomr_adv_dim1(idx) = V_flag(idx) * dXoutdX1
         dmomr_adv_di(idx)   = V_flag(idx) * dXoutdX2 + (z1-V_flag(idx)) * dXindX2
         dmomr_adv_dip1(idx) = V_flag(idx) * dXoutdX3 + (z1-V_flag(idx)) * dXindX3
         dmomr_adv_dip2(idx) =                          (z1-V_flag(idx)) * dXindX4


!-----------------------------------------------------------------------
!  Advektion des Strahlungsflusses
!-----------------------------------------------------------------------

         call ADV_SCHEME(XZ(MH,idx-1), XZ(MH,idx), XZ(MH,idx+1), XZ(MH,idx+2), &
                         Xout, dXoutdX1, dXoutdX2, dXoutdX3, Xin, dXindX2, dXindX3, dXindX4)

         H_adv(idx)          = V_flag(idx) * Xout     + (z1-V_flag(idx)) * Xin
         dH_adv_dim1(idx)    = V_flag(idx) * dXoutdX1
         dH_adv_di(idx)      = V_flag(idx) * dXoutdX2 + (z1-V_flag(idx)) * dXindX2
         dH_adv_dip1(idx)    = V_flag(idx) * dXoutdX3 + (z1-V_flag(idx)) * dXindX3
         dH_adv_dip2(idx)    =                          (z1-V_flag(idx)) * dXindX4

      end do
      !$OMP END PARALLEL DO



END SUBROUTINE ADVEC

END MODULE ADVECVAR

